function [R,F,Ab,LK] = resp2D(keff, in)
% new response function generator 2.23


dat   = in.dat;
abdat = in.abdat;
numg  = in.numg;
numm  = in.numm;
numel = in.numel;
order = in.order;

R    = zeros( in.numtypes, 4*numg*(order+1), 4*numg*(order+1) );
Rinc = zeros( numg*(order+1), numg*(order+1) );
Rfar = Rinc; Rlef=Rinc; Rrig=Rinc;
F    = zeros( in.numtypes, 4*numg*(order+1) );
Ab   = F;
tmpF = zeros( 1, numg*(order+1) );
tmpAb = tmpF;
LK = zeros( in.numtypes, 4,  numg*(order+1) );


% loop over all element types ---------------------------------------------
for el = 1:in.numtypes

mt  = in.mt(:,:,el);   

BCL = 0;
BCR = 0;
BCB = 2;
BCT = 0;      

src(:,:,1) = [ 0 
               0 ];            
src(:,:,2) = [ 0
               0 ];           

IBSL = [ 0 
         0 ]';
IBSR = [ 0
         0 ]';
IBSB = [ 1 
         0 ]';
IBST = [ 0 
         0 ]';

% loop over legendre order and group --------------------------------------
for RG = 1:numg     
for RO = 0:order

% -------------------------------------------------------------------------
[AB, AL, AC, AR, AT, S, dx, dy, dv, dc, siga, nusig, sct] = ...
    twoDcoefMGresp( dat, abdat, numg, numm, in.xcm, in.xfm, in.ycm, in.yfm, mt, src, 1, ...
	[BCL BCR BCB BCT], RO, RG );

N = sum(in.xfm); M = sum(in.yfm);
for i = 1:N+1
    for j = 1:M+1
        k = i+(j-1)*( N+1 );
        indx(k,1)=i;
        indx(k,2)=j;
    end
end

xint = zeros(N+1,1); 
for i = 2:N+1
    xint(i)=xint(i-1)+dx(i-1);
end
% pad end of AB,AL
AB = [AB' zeros(numg,N+1)]';
AL = [AL' zeros(numg,1)]';
% pad top of AR,AT
AT = [zeros(numg,N+1) AT']';
AR = [zeros(numg,1) AR']';
phi = zeros( (N+1)*(M+1), numg );
scsrc = zeros( (N+1)*(M+1), 1);
%==========================================================================
% set A1 and A2
% note to self -- will want to generalize this, maybe w/ a struct
A1 = spdiags([AB(:,1) AL(:,1) AC(:,1) AR(:,1) AT(:,1)], ...
    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
A2 = spdiags([AB(:,2) AL(:,2) AC(:,2) AR(:,2) AT(:,2)], ...
    [-N-1 -1 0 1 N+1],(N+1)*(M+1),(N+1)*(M+1));
% solve the first group flux
phieps = 1e-12;  mxit = 200;  it = 0; phierr = 1;

while (phierr > phieps && it < mxit)
phi0 = phi;
phi(:,1) = A1 \ ( S(:,1) +  1/keff * nusig(:,2).*phi(:,2) );
for g = 2:numg
    scsrc   = zeros( (N+1)*(M+1), 1);
    for gg = 1:g-1
        scsrc(:,1) = scsrc(:,1) + sct(:,g,gg).*phi(:,gg);
    end
    phi(:,g) = A2 \ ( S(:,g) + scsrc(:,1) );
end
phierr = norm( phi0 - phi );
it = it + 1;
end
if (it>=mxit)
    disp('exceeded mxit...')
    phierr
    disp('--')
end
%==========================================================================

% responses of various sorts, reponse due to  RG,  RO incident current
[CRL CRR CRB CRT LL LR LB LT] = cresp(phi,dx,dy,dc,mt,in.xfm,in.yfm,numg);
f1 = sum( phi(:,1).*nusig(:,1) );
f2 = sum( phi(:,2).*nusig(:,2) );
tmpF(1, (RO)*numg+RG) = f1+f2;
a1 = sum( phi(:,1).*siga(:,1) );
a2 = sum( phi(:,2).*siga(:,2) );
tmpAb(1, (RO)*numg+RG) = a1+a2;
% compute leakage
LK(el,1,(RO)*numg+RG) =  LB; %incident
LK(el,2,(RO)*numg+RG) =  LT; %far
LK(el,3,(RO)*numg+RG) =  LR; %to right of incident
LK(el,4,(RO)*numg+RG) =  LL; %to left of incident

% determine legendre coefficients -----------------------------------------
% Lp's should be put in a function and moved outside this loop
N = sum(in.xfm);    % max degree
k = 0:N;         % interval
maxord = order;  % maxorder to which we compute P's
P = zeros(N+1,N+1);
P(1,1:N+1) = 1; 
P(2,1:N+1) = 1 - (2*k)/N;
for j = 3:N+1 %otherwise N-1
    i = j-2;
    P(i+2,1:N+1) = ( (2*i+1)*(N-2*k).*P(i+1,:) - i*(N+i+1)*P(i,:) ) ./ ...
        ((i+1)*(N-i)) ;
end
for i = 0:maxord 
    a = (2*i+1);
    b = factorial(N+i+1)/factorial(N);
    c = factorial(N)/factorial(N-i);
    w(i+1)=a*c/b;
end
x = xint;

for G=1:numg
    % f def'n ASSUMES four-way symmetry
    f(1,:) = CRB(end:-1:1,G)';
    f(2,:) = CRT(:,G)';  
    f(3,:) = CRR(end:-1:1,G)';
    f(4,:) = CRL(:,G)';
 
    % now we want to expand f into fi
    for ord=1:maxord+1
        ff(:,ord) = w(ord) * f(:,:) * P(ord,:)';
    end
    side=1;
    Rinc( G:numg:end, (RO)*numg+RG ) = ff(1,:)';
    Rfar( G:numg:end, (RO)*numg+RG ) = ff(2,:)';
    Rrig( G:numg:end, (RO)*numg+RG ) = ff(3,:)';
    Rlef( G:numg:end, (RO)*numg+RG ) = ff(4,:)';
end

% end loop over legendre order and group ----------------------------------
end
end
% -------------------------------------------------------------------------
R(el,:,:) = [  Rinc  Rfar  Rrig  Rlef
               Rfar  Rinc  Rlef  Rrig
               Rlef  Rrig  Rinc  Rfar
               Rrig  Rlef  Rfar  Rinc ];               
F(el,:)   = [ tmpF tmpF tmpF tmpF];
Ab(el,:)  = [ tmpAb tmpAb tmpAb tmpAb];

% end loop over element types ---------------------------------------------
end

% % build total response matrix
% % this gets the info for each element type; now, need to place
% I = length(in.elements(:,1));
% J = length(in.elements(1,:));
% k = 0;
% % this should go in a function, or at the end of resp2D;
% for j = 1:J      % rows
%     for i = 1:I  % cols
%         if in.elements(i,j) > 0
%             % I am the kth element
%             k = k + 1;
%             type = in.elements(i,j);
%         else
%             % else I'm just void
%             continue
%         end
%         R( ((k-1)*lenR+1):(k*lenR),((k-1)*lenR+1):(k*lenR) ) = Rb(type,:,:);
%         F(1, ((k-1)*lenR+1):(k*lenR) ) = Fb(type,:);
%         A(1, ((k-1)*lenR+1):(k*lenR) ) = Ab(type,:);
%     end
% end



